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Abstract 

We study the critical behavior of the three-dimensional ± J Ising model [with a random-exchange 
probability P(J xy ) = p5(J xy — J) + (l—p)5(J xy + J)] at the transition line between the paramagnetic 
and ferromagnetic phase, which extends from p = 1 to a multicritical (Nishimori) point at p = 
Pn ~ 0.767. By a finite-size scaling analysis of Monte Carlo simulations at various values of p 
in the region pjy < p < 1, we provide strong numerical evidence that the critical behavior along 
the ferromagnetic transition line belongs to the same universality class as the three-dimensional 
randomly-dilute Ising model. We obtain the results v = 0.682(3) and n = 0.036(2) for the critical 
exponents, which are consistent with the estimates v = 0.683(2) and n = 0.036(1) at the transition 
of randomly-dilute Ising models. 

PACS numbers: 75.10.Nr, 75.40.Cx, 75.40.Mg, 64.60.Fr 
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I. INTRODUCTION 



The ± J Ising model has played an important role in the study of the effects of quenched 
random disorder and frustration on Ising systems. It is defined by the lattice Hamiltonian 

7~L±J = — ^ JxyCx^y, (1) 

where a x = ±1, the sum is over the nearest-neighbor sites of a simple cubic lattice, and the 
exchange interactions J xy are uncorrelated quenched random variables, taking values ±J 
with probability distribution 

P(Jxy) = p5(J X y ~J) + (l- P)S(J XV + J). (2) 

For p = 1 we recover the standard Ising model, while for p = 1/2 we obtain the usual 
bimodal Ising spin-glass model. 

The phase diagram of the three-dimensional (3D) ±J Ising model is sketched in Fig. 1. 
The high-temperature phase is paramagnetic for any p. The low-temperature phase depends 
on the value of p: it is ferromagnetic for small values of 1 — p, while it is a spin-glass 
phase with vanishing magnetization for sufficiently large values of 1 — p. The different 
phases are separated by transition lines, which meet at a multicritical point N located 
along the so-called Nishimori line. 1 ' 2,3 The spin-glass transition has been mostly studied at 
the symmetric point p = 1/2, see, e.g., Refs. 3,4 and references therein. The spin-glass 
transition line extends up to the Nishimori multicritical point, 2 located at 5,6 ' 7,8 p N 0.767. 
For larger values of p, the transition is ferromagnetic, up to p = 1 where one recovers the pure 
Ising model, and therefore a transition in the Ising universality class. At the ferromagnetic 
transition line, for p^ < p < 1, the critical behavior is expected to belong to a different 
universality class. 

An interesting hypothesis, which has already been put forward in Refs. 3,9, is that the 
ferromagnetic transition of the ±J Ising model belongs to the 3D randomly-dilute Ising 
(RDIs) universality class (see, e.g., Refs. 10,11 for reviews on randomly-dilute spin mod- 
els). A representative of the RDIs universality class is the randomly site-dilute Ising model 
(RSIM) defined by the lattice Hamiltonian 

H d = -J ^2 p x p y a x a y , (3) 

<xy> 
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FIG. 1: Sketch of the phase diagram of the 3D dbJ Ising model in the T-p plane. 

where p x are uncorrelated quenched random variables, which are equal to 0, 1 with proba- 
bility 

P( Px )=p6(p x -l) + (l-p)6( Px ). (4) 

Forp < 1 and above the percolation threshold of the spins (p pe rc ~ 0.3116081(13) on a cubic 
lattice 12 ), the RSIM under goes a continuous phase transition between a disordered and a 
ferromagnetic phase, whose nature is independent of p. This transition is definitely different 
from the usual Ising transition: for instance, the correlation-length critical exponent 13 ' 14 ' 15 ' 16 
v = 0.683(2) differs from the Ising value 17 ' 18 v = 0.63012(16). The RDIs universality class is 
expected to describe the ferromagnetic transition in generic diluted ferromagnetic systems. 
For instance, it has been verified that also the randomly bond-diluted Ising model (RBIM) 
belongs to the RDIs universality class. 13 ' 19 These results do not necessarily imply that also 
the ±J Ising model has an RDIs ferromagnetic transition line. Indeed, while the RSIM (3) 
has only ferromagnetic exchange interactions, the ± J Ising model is frustrated for any value 
of p < 1. Therefore, the ferromagnetic transition in the ±J Ising model belongs to the RDIs 
universality class only if frustration is irrelevant, a fact that is not obvious and should be 
carefully investigated. 

Reference 9 investigated the issue by means of a Monte Carlo (MC) renormalization- 
group (RG) study, claiming that the ± J Ising model belongs to the same RDIs universality 
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class as the RSIM and the RBIM. It should be noted however that the quoted estimate for 
the correlation-length exponent at the ferromagnetic transition, v = 0.658(9), is close to 
but not fully consistent with the RDIs value v = 0.683(2). 13 Another numerical MC work 20 
investigated the nonequilibrium relaxation dynamics of the ±J Ising model and showed an 
apparent nonuniversal dynamical critical behavior along the ferromagnetic transition line. 
These results are not conclusive and further investigation is called for to clarify this issue. 

In this paper we focus on the transition line of the 3D ±J Ising model between the 
paramagnetic and the ferromagnetic phase. We investigate the critical behavior by means 
of MC simulations at various values of p in the region p N < p < 1. Our finite-size scaling 
(FSS) analysis provides a strong evidence that the critical behavior of the 3D ±J Ising 
along the ferromagnetic line belongs to the 3D RDIs universality class. For example, we 
obtain v = 0.682(3) and rj = 0.036(2), which are in good agreement with the presently most 
accurate estimates 13 v = 0.683(2) and r\ = 0.036(1) for the 3D RDIs universality class. 

The paper is organized as follows. In Sec. II we summarize some FSS results which are 
needed for the analysis of the MC data, and describe our strategy to check whether the 
transition belongs to the RDIs universality class. In Sec. Ill we describe the MC simula- 
tions. In Sec. IV we report the results of the FSS analysis. Finally, in Sec. V we draw our 
conclusions. In App. A we report the definitions of the quantities we compute. 

II. STRATEGY OF THE FINITE-SIZE SCALING ANALYSIS 

In this work we check whether the ferromagnetic transition line in the 3D ± J Ising models 
belongs to the RDIs universality class. For this purpose, we present a FSS analysis of MC 
data for various values of p in the region 1 > p > p N « 0.767. We follow closely Ref. 13, 
which studied the ferromagnetic transition line in the 3D RSIM and RBIM and provided 
strong numerical evidence that these transitions belong to the same RDIs universality class. 
We refer to Ref. 13 for notations (a short summary is reported in App. A) and a detailed 
discussion of FSS in these disordered systems. 

According to the RG, in the case of periodic boundary conditions and for L — > oo, where 
L is the lattice size, a generic RG invariant quantity R at the critical temperature l/f3 c 
behaves as 

R(L, (3 = (3 C ) = R* (1 + c u L- u + c 12 L- 2u + ■■■ + c^L'^ + •••), (5) 
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where R* is the universal infinite- volume limit and uo and uj 2 are the leading and next- 
to-leading correction-to-scaling exponents. In RDIs systems scaling corrections play an 
important role, 16 ' 21 since u is quite small. Indeed we have tu = 0.33(3) and uo 2 = 0.82(8) 
in the 3D RDIs universality class. 13 These slowly-decaying scaling corrections make the 
accurate determination of the universal asymptotic behavior quite difficult. 

Instead of computing the various quantities at fixed Hamiltonian parameters, we keep 
a RG invariant quantity R fixed at a given value Rf 22 This means that, for each L, we 
determine the pseudocritical inverse temperature f3f(L) such that 

R(J3 = P f (L),L) = R f . (6) 

All interesting thermodynamic quantities are then computed at (5 — f3f(L). The pseudocrit- 
ical inverse temperature f3f(L) converges to (3 C as L — > 00. The value Rf can be specified 
at will, as long as Rf is taken between the high- and low-temperature fixed-point values 
of R. The choice Rf = R* (where R* is the critical-point value) improves the conver- 
gence of j3f to P c for L — * 00; indeed (3f — (3 C = 0{L~ 1 / V ) for generic values of Rf, while 
P f -fi c = 0(L~ l / v ~ u ) for Rf = R*. This FSS method has already been applied to the study 
of the critical behavior of N- vector spin models, 22 ' 23 and of randomly-dilute Ising models. 13 
As in Ref. 13, we perform a FSS analysis at fixed i?g = £/L = 0.5943, which is very close to 
the fixed-point value i?| = 0.5944(7) of R^ at j3 c . Given any RG invariant quantity R, such as 
the quartic cumulants U4 and U22, we consider its value at fixed R^, i.e., R(L) = R(L, f3f(L)). 
For L — > 00, R(L) behaves as R(L, f3 c ): 

R(L) =R*(1 + bnL'" + b 12 L- 2uJ + ■■■ + b 21 L~^ + •••), (7) 

where the coefficients depend on the Hamiltonian. The derivative R! with respect to (3 
of a generic RG invariant quantity R behaves as 

R'(L) = aL l/v (1 + anL-" + a 12 L~ 2uJ + ■■■ + a 2l L~^ + •••). (8) 

Finally, the FSS of the magnetic susceptibility x is given by 13 

X (L) = X (L, P = Pf(L)) = el 2 ^ (l + e 11 L~ u + e l2 L^ + ■■■ + e 21 L~^ + •••)+ e b (9) 

where represents the background contribution. 

A standard RG analysis, see, e.g., Ref. 13, shows that the amplitudes of the 0(L~ kul ) 
scaling corrections are proportional to (with a universal coefficient), where M3 is the 
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leading irrelevant scaling field with RG dimension y 3 = —to. Hamiltonians such that = 0— 
we call them improved Hamiltonians — have a faster approach to the universal asymptotic 
behavior, because the 0(L~ ku ) scaling corrections vanish: b lk = a lk = e lk = in Eqs. (7), 
(8), and (9). In this case the leading scaling corrections are proportional to -u 4 L~ W2 , where w 4 
is the next-to-leading irrelevant scaling field and 2/4 = — lo<i is its RG dimension. In Ref. 13 
it was shown that the RSIM for p = p* = 0.800(5) and the RBIM for p = p* = 0.56(2) 
are improved. Since scaling fields are analytic functions of the Hamiltonian parameters, M3 
must be proportional to p — p* close to p = p*, i.e. w 3 ~ — p*). Therefore, since the 
coefficients bi k , a lfc , and e lfc that appear in Eqs. (7), (8), and (9) are proportional to m|, we 
have 

h k ,a lk ,e lk ~ (p -p*) k . (10) 

Beside the quantities defined in App. A, we also consider observables — in analogy with 
the previous terminology, we call them improved quantities — characterized by the fact that 
the leading scaling correction proportional to L~ w (approximately) vanishes in any model 
belonging to the RDIs universality class. 13 We consider the combination of quartic cumulants 

Uun = UA + 1-3C/22, (11) 

and improved estimators of the critical exponent v defined as 

^,im = R'fl U' 4;im = U'jj™ (12) 

(Ud is defined in App. A). In Ref. 13 we showed that, if the transition belongs to the RDIs 
universality class, the leading scaling correction proportional to of these improved ob- 
servables is suppressed. More precisely, we showed that the universal ratio of the amplitudes 
of the leading scaling correction in U im and U4 satisfies 

\bn,ujb n ,u 4 \ < ^, (13) 

while the one for the quantities R'^ im and R'^ is bounded by 

1 

Wi^Ja u ^\<-. (14) 

The remaining scaling corrections are of order L~ 2uJ and L~ U2 . These improved observables 
are particular useful to check whether the transition in a given system belongs to the 3D 
RDIs universality class. 
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FIG. 2: Exponential autocorrelation time r of the magnetic susceptibility for a mixture of 
Metropolis and cluster updates as discussed in the text, at p = 0.87. The dotted line shows 
the result of a fit to r = cL z : this fit gives z ~ 1.6. 

To summarize: in order to check whether the ferromagnetic transition of the 3D ±J 
Ising model belongs to the RDIs universality class, we perform a FSS analysis at fixed 
i?g = 0.5943, and check if the results for the critical exponents and other universal quantities 
are consistent with those obtained for the RDIs universality class, which is characterized 
by 13 critical exponents v = 0.683(2) and rj = 0.036(1), by the leading and next-to-leading 
scaling-correction exponents u = 0.33(3) and 0J2 = 0.82(8) and by the universal infinite- 
volume values of the quartic cumulants LV" 2 * 2 = 0.148(1), U* m = 1.840(4), and XJ* d = 1.500(1). 
Notice that the fact that we fix = 0.5943 does not introduce any bias in our FSS analysis. 

III. MONTE CARLO SIMULATIONS 

We performed MC simulations of Hamiltonian (1) with J = 1 for p = 
0.94, 0.90, 0.883, 0.87, 0.83, 0.80, close to the critical temperature on cubic lattices of size 
L 3 with periodic boundary conditions, for a large range of lattice sizes: from L = 8 to 
L = 80 for p = 0.883, 0.87, to L = 64 for p = 0.94, 0.90, 0.83, and to L = 48 for p = 0.80. We 
chose values of p not too close to p = 1: indeed, as p — > 1 we expect crossover effects due to 
the presence of the Ising transition for p = 1 and, therefore, that the asymptotic behavior 
sets in only for large values of L. We return to this point later. 



We used a Metropolis algorithm and multispin coding. 24 In the simulation n^ t systems 
evolve in parallel, where nbit = 32 or nut = 64 depending on the computer that is used. 
For each of these n hit systems we use a different set of couplings J xy . This allows us to 
perform 64 parallel simulations on a 64-bit machine, and therefore to gain a large factor 
in the efficiency of the MC simulations. We used high-quality random-number generators, 
such as the RANLUX 25 or the twister 26 generators. 27 Using the twister random-number 
generator, we need about 1.2 x 10~ 9 seconds for one Metropolis update of a single spin on 
an Opteron processor running at 2 GHz. Our simulations took approximately 3 CPU years 
on an Opteron (2 GHz) processor. 

It is worth mentioning that cluster algorithms, such as the Swendsen-Wang cluster 28 
and the Wolff single-cluster 29 algorithm, show significant slowing down in the ±J Ising 
model. At the earlier stage of this work we performed some simulations of the ±J Ising 
model at p — 0.87 using the algorithm used in Ref. 13 to simulate the RSIM and the 
RBIM. There we used a combination of Metropolis, Swendsen-Wang cluster, 28 and Wolff 
single-cluster 29 updates. More precisely, each updating step consisted of 1 Swendsen-Wang 
update, 1 Metropolis update, and L single-cluster updates. In all cases the exponential 
autocorrelation times r of the magnetic susceptibility was small: r < 1 in units of the 
above updating step, even for the largest lattice sizes considered, i.e. L = 192. In the ±J 
Ising model at p — 0.87 autocorrelation times are much larger. In Fig. 2 we plot estimates 
of r as obtained from the magnetic susceptibility. They show a clear evidence of critical 
slowing down: r ~ L z with z ~ 1.6. Such a value of z should be compared with the 
dynamic exponent of Swendsen-Wang and Wolff cluster algorithms in the RSIM, which is 
much smaller: 30 z < 0.5. These results show that cluster algorithms behave differently in 
the ±J Ising model, likely due to frustration. They suggest that frustration is relevant for 
the cluster dynamics. 

Taking also into account the computer time required by the cluster algorithms, we then 
turned to a multispin Metropolis algorithm. This turns out to be much more effective at the 
lattice sizes considered, although it has a larger dynamic exponent z > 2, see, e.g., Ref. 30 
and references therein. We also mention that the autocorrelation time significantly increases 
with decreasing p (keeping L fixed). For example, for L = 48 it increases by approximately 
a factor of 10 from p = 0.90 to p — 0.80. This represents a major limitation to perform 
simulations for large lattices close to the multicritical point. 
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For each lattice size we considered N s disorder samples, with N s decreasing with in- 
creasing L, from N s > 10 6 for L = 8 to N s > 2 x 10 4 for the largest lattices. For each 
disorder sample, we collected a few hundred independent measurements at equilibrium. The 
averages over disorder are affected by a bias due to the finite number of measures at fixed 
disorder. 13,31 A bias correction is required whenever one considers the disorder average of 
combinations of thermal averages. We used the formulas reported in App. B of Ref. 13. 
Errors were computed from the sample-to-sample fluctuations and were determined by using 
the jackknife method. 27 

Our FSS analysis is performed at fixed = £/L. In order to determine expectation 
values at fixed R%, one needs the values of the observables as a function of j3 in some 
neighborhood of the inverse temperature (5 Tun used in the simulation. In Ref. 13 we used the 
reweighting method for this purpose. This requires that the observables and, in particular, 
the values of the energy are stored at each measurement. For the huge statistics like those we 
have for the smaller values of L, this becomes unpractical. Therefore, we used here a second- 
order Taylor expansion, determining 0((3, L) from 0(/3 mn , L) + a {(3 — /3 mn ) + b {f3 — /3 mn ) 2 . 
The coefficients a Q and bo are obtained from appropriate expectation values as in Ref. 23. 
Since their computation involves disorder averages of products of thermal averages, we have 
implemented in all cases an exact bias correction, using the formulas of Ref. 13. Derivatives 
with respect to (3 are then obtained as ; (f3,L) = ao + 2&o(/3 — A-un)- Of course, this 
method requires |/3 run — f3f\ to be sufficiently small. We have carefully checked the results 
by performing, for each L and p, runs at different values of (3. 

The MC estimates of the quantities introduced in Sec. II and in App. A at fixed R^ = 
i/L = 0.5943 are available on request. 

IV. FINITE-SIZE SCALING ANALYSIS 

In this section we present the results of our FSS analysis of the MC data at fixed R^ = 
0.5943. 
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FIG. 3: MC estimates of U22 versus L u with u = 0.33 for different values of p. The dotted lines 
show results of fits to Co + c\L~ £l + C2L~ £2 , fixing cq = 0.148, £\ = 0.33, and £2 = 0.82. In the 
RDIs universality class C/ 2 * 2 = 0.148(1). 13 

A. Renormalization-group invariant quantities 

In Fig. 3 we show the MC estimates of U 22 versus L _£J with uo = 0.33(3), which is the 
leading scaling exponent of the RDIs universality class. The data vary significantly with p 
and L. This p and L dependence is always consistent with the existence of the expected 
next-to-leading scaling corrections, i.e. with a behavior of the form 

U22 = U22 + ciL- £l + c 2 L~ £ \ (15) 

where C7| 2 , £1 arid e 2 are fixed to the RDIs values: 13 = 0.148, E\ = 0.33 and e 2 = 
0.66, 0.82. The fits corresponding to 62 = 0.82 are shown in Fig. 3. Note that in most of 
the cases it is crucial to include a next-to-leading correction. Only for p = 0.90 the data are 
well fitted by taking only the leading scaling correction. 

An unbiased estimate of 00 can be obtained from the difference of data at different values 
of p, i.e. by considering 

U 22 ( Pl ;L)-U 22 (p 2 ;L)ncL-". (16) 
Linear fits of the logarithm of these differences give results in reasonable agreement with 
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FIG. 4: MC estimates of [/i m versus L 0,82 . The filled square on the vertical axis corresponds to 
the RDIs estimate 13 L^* m = 1.840(4). 

the RDIs estimate u = 0.33(3), especially when only data corresponding to L > L min = 24 
are used. For L min = 24 [L min = 32], we obtain to = 0.27(2) [a; = 0.27(3)] from the data at 
pi = 0.83 and p 2 = 0.90, u = 0.19(5) [u = 0.31(9)] from those at Pl = 0.883 and p 2 = 0.90, 
and u = 0.29(3) [u = 0.25(4)] from the results at Pl = 0.83 and p 2 = 0.883. We also fitted 
the difference U22 — 0.148 at p — 0.90 to cL~ e (for this value of p next-to-leading corrections 
are apparently very small, see Fig. 3). We obtain u = 0.35(3) [uj = 0.39(6)] for L min = 24 

[-^min = 32]. 

The results of the above-reported fits of U22 show that the leading scaling corrections 
proportional to vanish for p rs 0.883. Note that, close to p*, the relevant next-to- 
leading scaling corrections should be those proportional to L - " 2 with 002 ~ 0.82. Indeed, 



according to Eq. (10), the coefficient of those proportional to L~ 2uJ is of order (p — p*) 2 , 
i.e. bi2 ~ bu^p — p*). Therefore, 612 is small if 612 = 0(1) (we checked this numerically). 
This applies to the FSS at p — 0.87 and 0.90, where the L~ 2uJ corrections can be neglected, 
although in these two cases we cannot neglect the leading correction whose coefficient is 
proportional to p — p*. An analysis of the leading scaling corrections at p = 0.87, 0.883, 0.90, 
assuming the RDIs values TJ£ 2 = 0.148(1) and uj = 0.33(3) (we perform combined fits to (15) 
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FIG. 5: Estimates of U* m as obtained by fits to U* m + cL~ £ , versus the minimum lattice size L m , n 
allowed in the fits. Some data are slightly shifted along the x-axis to make them visible. The 
dotted lines correspond to the RDIs estimate 13 U* = 1.840(4). 

with ei = to) gives the estimate 

p* = 0.883(3), (17) 

which is approximately in the middle of the ferromagnetic line, i.e. 1 — p* ~ (1 — p^)/2. 
We performed a similar analysis for C/ d) obtaining a consistent estimate of p*. Thus, the ±J 
Ising model for p = 0.883 is approximately improved. Therefore, at p — 0.883, fits of the 
data assuming 0(L - ^ 2 ) leading scaling corrections should provide reliable results. 

As discussed in Sec. II, a useful quantity to perform stringent checks of universality within 
the RDIs universality class is the combination U- im of quartic cumulants reported in Eq. (11). 
For this quantity the scaling corrections proportional to L~ w are small, cf. Eq. (13), and 
thus the dominant corrections should behave as L~ 2uJ , with 2u ~ 0.66. As already discussed, 
for values of p close to p*, such as p = 0.87, 0.883, 0.90, also the L~ 2uJ term is expected to be 
small and thus the dominant corrections should scale as L~ W2 with 002 ~ 0.82. In Fig. 4 we 
show the MC results for U irn for various values of p. Fig. 5 shows results of fits to 

U* + cL- £ , (18) 
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FIG. 6: Estimates of as obtained by fits of Ud at p = 0.883, to + cL 6 . Some data are 
slightly shifted along the x-axis to make them visible. The dotted lines correspond to the RDIs 
estimate 13 E/J = 1.500(1). 

with e = 0.66,0.82. We obtain U* m = 1.840(3) [3], 1.842(2) [1], 1.845(2) [3] respectively for 
p = 0.90, 0.883, 0.87, fixing e = uj 2 = 0.82(8) (the error in brackets is related to the uncer- 
tainty of lo 2 ) and using data with L > 32; moreover we obtain U* m = 1.847(3) [2], 1.840(7) [1] 
respectively for p = 0.83, 0.80, fixing e = 0.66(6) and using data with L > 24. For all values 
of p the results are in good agreement with RDIs estimate 13 U* m = 1.840(4). They provide 
strong support to a RDIs critical behavior along the ferromagnetic line. 

A further stringent check of universality comes from the analysis of the data Ud at 
p = 0.883, because the data of Ud are very precise due to a cancellation of the statisti- 
cal fluctuations. 13 Since the model is improved, the L~ kuj scaling corrections are negligible 
and the large-L behavior is approached with corrections of order L'^ 2 , uj 2 = 0.82(8). We 
thus fit the data to Eq. (18) with s = 0.74,0.82,0.90. In Fig. 6 we show the results. We 
obtain U* d = 1.5001(1) [15], 1.5004(2) [9], 1.5003(10) [6] (the error in brackets is related to the 
uncertainty of ui 2 = 0.82(8)) for L min = 12, 24, 48 respectively. Moreover, by fitting the data 
to £/ d * + ciL- ei +c 2 L~ e2 with^ = 0.33 and e 2 = 0.82, we obtain U* = 1.5006(7), 1.500(3) for 
L m in — 12, 24 respectively. These results are in perfect agreement with the RDIs estimate 13 
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FIG. 7: Estimates of the critical exponent u, as obtained by fits of R'g, U 4 , R'^ im , and U' 4 im at 
p = 0.883. L m [ n is the minimum lattice size allowed in the fits. Some data are slightly shifted along 
the x-axis to make them visible. The dotted lines correspond to the estimate v = 0.682(3). 

U*[ = 1.500(1). Such an agreement is also confirmed by the analysis of the data of U22, for 
example a fit to Eq. (18) with e = uj 2 = 0.82(8) gives C/* 2 = 0.1486(8) [3] for L min = 32, to 
be compared with the RDIs estimate 13 [7| 2 = 0.148(1). 

We have not shown results for values of p too close to 1, for p > 0.90 say, because they 
are affected by crossover effects due to presence of the Ising transition for p = 1, as it also 
occurs in randomly dilute Ising models. 16,21,32 For instance, for p = 0.94 the data are not 
compatible with a behavior of the form (15) with U% 2 fixed to the RDIs value. Our data 
that correspond to lattice sizes L < 64 apparently converge to a smaller value, consistently 
with the expected crossover from pure to random behavior (in pure systems C/| 2 = 0). The 
same quantitative differences are observed in the RSIM and in the RBIM close to the Ising 
transition. This suggests that in FSS analyses up to L w 100 the asymptotic RDIs behavior 
can only be observed for p < 0.94. 
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B. Critical exponents 



The correlation-length exponent v can be estimated by fitting the derivative of and U4 
to the expression (8). Accurate estimates are only obtained for improved Hamiltonians. For 
generic models, as shown in Ref. 13, good estimates are only obtained by using improved 
estimators, such as those reported in Eq. (12). 

We analyze the data at p — 0.883, which is a very good approximation of the improved 
value p* = 0.883(3). In Fig. 7 we report several results for the critical exponent u, obtained 
by analyzing R'^, U' A , and their improved versions R'^ im and U' A ira . We show results of fits of 
their logarithms to 

hnL + a + bL~ £ , (19) 

fixing e to several values. Since the Hamiltonian is approximately improved, scaling cor- 
rections are expected to decrease as L^ 2 with 002 = 0.82(8). Since p = 0.883 is only 
approximately equal to p*, one may be worried of the residual leading scaling corrections 
that are small but do not vanish exactly. Improved estimators should provide the most 
reliable results since the leading scaling corrections are additionally suppressed. 

As can be seen in Fig. 7, the results obtained by using R'^ and U' 4 and e = 0.82 are 
perfectly consistent with those obtained from their improved versions. This confirms that 
the Hamiltonian is improved. Fits of R'^ to (19) with e = 0.33 do not provide stable results. 
The results approach the values obtained in the other fits only when increasing the minimum 
size L min allowed in the fit. This is expected, since the L~ w corrections should be negligible 
with respect to the L^ 2 ones. In conclusion, our final estimate of the correlation-length 
exponent is 

v = 0.682(3), (20) 

which includes all results (with their errors) of the fits of R'^, R'^ im , U' 4 im to Eq. (19) 
with e = 0.82(8) and L min = 16, 24. Estimate (20) is in perfect agreement with the most 
precise RDIs estimate v = 0.683(2). 

Estimate (20) is also confirmed by the analysis of the data at the other values of p. Fig. 8 
shows results obtained by fitting the logarithm of R'^ im to the function (19) for other values 
of p. They are definitely consistent with the result obtained at p — 0.883. Results for 
p = 0.80 are not shown because the available data are not sufficient to get reliable results. 
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FIG. 8: Results from v obtained by fitting BL.^ to aL l / v (l + bL' £ ). Some data are slightly shifted 
along the x-axis to make them visible. The dotted lines correspond to the estimate obtained at 
p = 0.883, i.e. v = 0.682(3). 



In order to estimate the critical exponent 77, we analyze the FSS of the magnetic suscep- 
tibility x-i °f- Eq. (9). We fit it to aL 2 ~ v + b (where b represents a constant background 
term), to aL 2_,? (l + cL~ £ ), and to aL 2 ~ v (l + cL~ e ) + b (more precisely, we fit ln% to the 
logarithm of the previous expressions). The results at p — 0.883 are shown in Fig. 9, versus 



the minimum size L min allowed in the fits. We obtain the estimate 



77 = 0.036(2), 



(21) 



which includes all results obtained for L min > 16. This estimate agrees with the most precise 
RDIs estimate 77 = 0.036(1). Fig. 10 shows results for the other values of p. Again, they are 
in good agreement. 



C. The critical temperature 

The critical temperature can be estimated by extrapolating the estimates of (3f at = 
0.5943, cf. Eq. (6). Since we have chosen = 0.5943 w R% = 0.5944(7), 13 we expect 
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FIG. 9: Estimates of the critical exponent rj, obtained by fitting lnx at p = 0.883, to a + (2 — 
rj) In L + 6L^ 2 (denoted by +cost) , a + (2 - ry) In L + cti L~ £ , and to a + (2 - 77) In L + m L~ e + 6L^- 2 . 
Some data are slightly shifted along the x-axis to make them visible. The dotted lines correspond 
to the final estimate rj = 0.036(2). 

in general that (5f — f3 c — 0{L^ 1 l v ^ w '). For p = 0.883, since the model is approximately 
improved, the leading scaling corrections are related to the next-to-leading exponent u>2- 
Thus, in this case (3f — (3 C = 0{L~ x l v ~ w ' 1 '). This behavior is nicely observed in Fig. 11, which 
shows p f (L) at p = 0.883 vs L~ 1 / U ~ UJ ' 2 with 1/u + oo 2 ~ 2.28. A fit to (3 C + al" 1 /"^ 2 gives 
P c = 0.300611(1). 

For the other values of p we expect (3 f - (3 C = 0(L~ 1/! ^) with l/u + uo w 1.79. Linear 
fits of flf(L) (for L > L min with L min sufficiently large to give an acceptable \ 2 ) gi ye the 
estimates (3 C = 0.25544(2) for p = 0.94, /3 C = 0.285285(5) for p = 0.90, (3 C = 0.313748(1) for 
p = 0.87, p c = 0.365459(5) for p = 0.83, (3 C = 0.42501(3) for p = 0.80. We finally recall 
that 18 p c = 0.22165452(8) for p = 1 (the standard Ising model), and that 7 (3 C = 0.5967(11) 
at the multicritical Nishimori point at p^ = 0.7673(3). In Fig. 12 we plot the available 
estimates of the critical temperature T c = 1//3 C in the region 1 > p > p^. 

The estimates of T c shown in Fig. 12 hint at a smooth linear behavior for small values of 
w = 1 — p, close to the Ising point at w — 0. This can be explained by some considerations 
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FIG. 10: Estimates of the critical exponent 77, obtained by fitting x for various values of p. See the 
caption of Fig. 9 for an explanation of the fits. Some data are slightly shifted along the a>axis to 
make them visible. The dotted lines correspond to the result 77 = 0.036(2) obtained at p = 0.883. 

on the multicritical behavior around the Ising point at w = 0. The Ising critical behavior 
at w = is unstable against the RG perturbation induced by quenched disorder at w > 0, 33 
which leads to the RDIs critical behavior. Indeed such a perturbation has a positive RG 
dimension y w at the Ising fixed point: 10,34 y w = ais/^is = 2/Vi s — 3 where cvi s and v\ s are the 
Ising specific-heat and correlation-length critical exponents, and therefore 17 y w = 0.1740(8). 
Thus, in the absence of an external magnetic field, beside the scaling field u t related to the 
temperature, there is another relevant scaling field u w associated with the quenched disorder 
parameter w = 1 — p. General RG scaling arguments 21,35 show that the singular part of the 
free energy for w — > behaves as 

f si n g ~u t - ais F(X), X = u w u t \ (22) 

where <fi = y m v\ s = a\ s = 0.1096(5) is the crossover exponent, and F(X) is a crossover 
scaling function which is universal (apart from normalizations). The scaling fields Ut and 
u w depend on the parameters of the model. In general, we expect 

u t = t + kw, (23) 
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FIG. 11: Estimates of f3 f {L) at p = 0.883 versus ZrO+wa) for l/u + uj 2 ~ 2.28. The dashed line 
corresponds to a linear fit of the data for L > 12. 

where i = T/Tj s — 1, Ti s is the critical temperature of the Ising model, and k is a constant. 
No such mixing between t and u> occurs in u w , since w w vanishes for w = 0. Hence, we can 
take u w = w. The system has a critical transition for w > at T c (w). Since the singular 
part of the free energy close to a critical point behaves as (T — T c ) 2 ~ a (a = —0.049(6) is the 
specific-heat exponent of the RDIs universality class), we must have F(X C ) = 0, where X c 
is the value of X obtained by setting T = T c (w) (see, e.g., Ref. 36 and references therein). 
Hence, we obtain 

w [T c (w)/T ls - 1 + hv]~* = X c , (24) 

and therefore 

T c {w)/T ls - 1 = (w/X c )^ -kw + ---, (25) 

where the dots indicate higher-order terms. This expression provides the w dependence of 
the critical temperature for w small. Note that the nonanalytic term in Eq. (25) is suppressed 
with respect to the analytic ones, because 1/0 ~ 9.1. Thus, T c {w) ~ Ti s (l — kw + 0(w 2 )). 
Since T c (w) < Ti s , we can also infer that k > 0. From the results for T c (w) we estimate 
k pa 2.2 for the ± J Ising model. 
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FIG. 12: The critical temperature T c = 1/ f3 c vs 1 — p. 
V. CONCLUSIONS 



In this paper we have studied the critical behavior of the 3D ±J Ising model at the 
transition line between the paramagnetic and the ferromagnetic phase, which extends from 
p = 1 to a multicritical (Nishimori) point at p = pn ~ 0.767. We presented a FSS analysis of 
MC simulations at various values of p in the region p^ < p < 1. The results for the critical 
exponents and other universal quantities are consistent with those of the RDIs universality 
class. For example, we obtained v = 0.682(3) and rj = 0.036(2), which are in good agreement 
with the presently most accurate estimates 13 v = 0.683(2) and rj = 0.036(1) for the 3D RDIs 
universality class. Therefore, our FSS analysis provides a strong evidence that the critical 
behavior of the 3D ± J Ising along the ferromagnetic line belongs to the 3D RDIs universality 
class. 

We also note that the random-exchange interaction in the ±J Ising model gives rise 
to frustration, while the RDIs universality class describes transitions in generic diluted 
Ising systems with ferromagnetic exchange interactions. This implies that frustration is 
irrelevant at the ferromagnetic transition line of the 3D ±J Ising model. Moreover, the 
observed scaling corrections are consistent with the RDIs leading and next-to-leading scaling 
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correction exponents u> = 0.33(3) and 002 = 0.82(8). This indicates that frustration does not 
introduce new irrelevant perturbations at the RDIs fixed point with RG dimension y/ > — 1. 

APPENDIX A: NOTATIONS 

We define the two-point correlation function 



G(x) = (a a x ), (Al) 

where the overline indicates the quenched average over the J xy probability distribution. 
Then, we define the corresponding susceptibility x = ^2 X G( X ) anc ^ the correlation length £ 



„ 2 _ G(0) - G(q min ) 



(A2) 



?min^'(?min) 

where q min = (2n/L, 0,0), g = 2 sing/2, and G(q) is the Fourier transform of G(x). We also 
consider quantities that are invariant under RG transformations in the critical limit. Beside 
the ratio 

Ik = (A3) 
we consider the quartic cumulants U4, U22 and Ua defined by 

U, = (A4) 

U22 = —=5 , 

1^2 

U d = U A - U22, 

where 

X 

We also define corresponding quantities C7 4 , £7 22 , and C/ d at fixed = 0.5943. Finally, we 
consider the derivative R'^ of i?^, and U' 4 of U 4 , with respect to f3 = 1/T, which allow one to 
determine the critical exponent v. 
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